First-principles methodology for quantum transport in multiterminal junctions 
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We present a generalized approach for computing electron conductance and I-V characteristics 
in multiterminal junctions from first-principles. Within the framework of Keldysh theory, elec- 
tron transmission is evaluated employing an O(N) method for electronic-structure calculations. The 
nonequilibrium Green function for the nonequilibrium electron density of the multiterminal junction 
is computed self-consistently by solving Poisson equation after applying a realistic bias. We illus- 
trate the suitability of the method on two examples of four-terminal systems, a radialene molecule 
connected to carbon chains and two crossed carbon chains brought together closer and closer. We 
describe charge density, potential profile, and transmission of electrons between any two terminals. 
Finally, we discuss the applicability of this technique to study complex electronic devices. 

PACS numbers: 72.10.-d, 85.65.+h, 73.63.-b, 85.35.-p 
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I. INTRODUCTION 

Electron transport through molecular-scale devices has 
become a very exciting research area for both experimen- 
talists and theorists. The main reason for this interest 
originates from the possibility of extreme miniaturization 
in electronic devices. In recent years, hundreds of papers 
have been published to establish the connection between 
the microscopic characteristics of an electronic system, 
such as the atomic configuration and the electronic struc- 
ture, and transport properties such as electrical current 
and conductance. These results have helped to improve 
the understanding of the I-V characteristics of nanojunc- 
tions. However, the interpretation of the I-V curves in 
terms of the geometry of the junction remains largely a 
fundamental challenge for molecular electronics. In fact, 
it has not yet been possible to establish a general theo- 
retical model that can reliably deal with any molecular 
junction of arbitrary geometry. Owing to the complex- 
ity of the system, these studies arc strongly dependent 
on the existence of reliable theoretical treatments based 
on first-principles approaches. In some instances, even 
conventional approaches based on density-functional the- 
ory (DFT) are expected to break down, especially in the 
case of weak coupling^. Nevertheless, DFT is expected 
to provide reliable results in a large number of cases. To 
the best of our knowledge, all the existing approaches 
to date, based on ab initio calculations, can deal with 
systems limited to two terminals only. It is therefore 
of widespread interest to develop robust computational 
schemes that can routinely and reliably account for the 
transport mechanism in multiterminal molecular devices. 

In his seminal work, Biittiker— developed a conduc- 
tance formula for a four-terminal system. However, that 
formula was not explicitly implemented in the frame- 
work of first-principles based calculations. Moreover, 
since the idea of that paper was to propose a reliable 
method for voltage difference measurements, Biittiker as- 
sumed that only two of the leads can carry current to 



and from the sample and the two others only measure 
the voltage. A similar non-atomistic approach based 
on tight-binding approximation has been presented by 
Baranger et al£. Within the framework of a Luttingcr 
liquid theory the four-terminal resistance of an inter- 
acting quantum wire was studied by Arrachea et alA. 
Recently, Jayasckcra et al. proposed a four-terminal 
approach for magneto-transport properties based on R- 
matrix theory^. However the approach is only applica- 
ble to two-dimensional devices, it is formulated in the 
framework of semi-empirical tight-binding, and does not 
include a self-consistent treatment of finite applied po- 
tential. Finally, a mesoscopic treatment for phonon- 
assistcd current through multiterminal conductors was 
formulated by Rychkov et al&. While these multiter- 
minal approaches address important issues, they neither 
treat the system in an ab initio fashion, including atom- 
istic details, nor do they account for the self-consistent 
(SC) rearrangement of electrons as the bias and the cur- 
rent increase. The importance of the self-consistency had 
been demonstrated in our previous paper—, showing that 
negative differential resistance in the I-V characteristic 
can only be quantitatively studied when self-consistency 
is included. 

In this paper, we present a generalized approach for 
computing conductance and I-V characteristics in mul- 
titerminal junctions, based on density-functional the- 
ory. In order to take into account the difference in 
electro-chemical potentials in different leads, we use non- 
equilibrium Keldysh formalism. The electronic trans- 
port is formulated in the basis of an O(N) method for 
electronic-structure calculations, which is an ab initio 
pscudopotcntial density functional approach using a lin- 
ear combination of numerical atomic orbitals (LCAO) ba- 
sis that are optimized for the problem in hand^. We apply 
external bias voltage through any lead in a realistic way 
and self-consistently compute the nonequilibrium Green 
function for the nonequilibrium electron density of the 
multiterminal junction by solving Poisson equation. One 
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FIG. 1: (Color online) Schematic diagram of a multitermi- 
nal junction. The barrier region is divided into n-number of 
blocks Ci, C2, C3, ■ • ■ , C n . The leads Li, L2 are connected to 
blocks Ci, C n , respectively, whereas the remaining leads L3, L4 
are connected to block C3. 



of the main advantages of our scheme is that we can apply 
bias through any number of leads and at the same time 
compute current between any two leads. The method is 
illustrated on two prototypical four-terminal systems (i) 
a radialene molecule connected to carbon chains, and (ii) 
two crossed carbon chains brought together closer and 
closer. We discuss the charge density, potential profile, 
and transmission of electrons between any two terminals. 
We also evaluate the current flowing between them. The 
algorithmic approach has been implemented on massively 
parallel computer architectures, and can therefore be ap- 
plied to systems of realistic sizes. 

The paper is organized as follows. The basic the- 
ory is outlined in section [ill sketching a noncquilibrium 
Green function formulation of the conductance calcula- 
tions. The applications are discussed in Section HTT1 The 
scheme is applied to two simple systems, addressing in 
particular the symmetry in transmission and effect of the 
different bias voltages on the I-V curves. The paper con- 
cludes with a summary in Section [IVI 



II. THEORY 

In this section, we describe the theoretical aspects of 
conductance calculations. The general system setup, the 
coupling to the leads, the equilibrium and nonequilibrium 
density matrices, the implementation of the bias voltage 
through any number of leads, and the conductance for- 
mula are presented. 



A. System setup 

We consider the prototypical multiterminal system 
sketched in Fig. [TJ Two or more semi-infinite leads 
Li,L2,L3,... are coupled to a central barrier region 
C with thermal reservoirs that are maintained at the 
electro-chemical potentials fi%, fj,2, A*3) ' "• Within our 
approach, region C can be treated as n subregions as 
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FIG. 2: (Color online) (Left) Schematic of a four-terminal sys- 
tem. The leads L2, L3 and L4 are connected to the molecu- 
lar barrier M via the subregions Ci, C2, C3, and C4 respectively. 
The subregions C3, C4 and the molecular region M are consid- 
ered together as a single region C shown by the blue-dotted 
box. The black-dotted box shows the extended-scattering re- 
gion S. (Right) Tri-diagonal matrix representation of this 
system. Any matrix-element in the lower off-diagonal-blocks 
is the complex conjugate of the corresponding element in the 
upper off-diagonal-blocks and therefore one can avoid stor- 
ing the lower off-diagonal-blocks. The shaded blocks of the 
matrix are directly connected to the leads and hence these 
contain the same matrix-elements as in the respective leads. 



Ci, C2, C3, . . . , C n and, in principle, we may connect m- 
number of leads to them. Note that an important hy- 
pothesis of the approach (also implicit in all two-terminal 
approaches based on Green function) is that there are no 
direct interactions between the leads and they only in- 
teract via the barrier region. Consequently the overlap 
integrals between orbitals on atoms situated in different 
leads take place via the barrier region only. 



For reasons of simplicity, we will be discussing a four- 
terminal system as an example. However, the method 
can be readily generalized to any number of electrodes. 
Suppose the leads l_i, L2, L3 and L4 are connected to the 
molecular barrier M via the subregions Ci, C2, C3 and C 4 
respectively. At the beginning of the calculation, each 
subregion has the same potential and charge distribution 
as its respective connected lead. 



In order to study the transport properties, in principle, 
we need to invert an infinite Hamiltonian of the infinite 
system which includes all parts of the semi-infinite leads. 
However, the electrons injected from the reservoirs move 
ballistically through the leads and all scattering events 
only occur around the barrier region (molecular region), 
called extended-scattering region S [see Fig. |2] (left)]. The 
potential is modified only within a finite region of the 
leads, the one being connected to the barrier region. It is 
therefore sufficient to consider a finite Hamiltonian con- 
taining all subregions Ci, C2, C3, C4 and the barrier region 
M. The Hamiltonian matrix of this system is of finite 
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(1) 

where i?c» j -Hm are the Hamiltonian matrices in the ith- 
lead and the molecular barrier respectively, and Vqm is 
the interaction between the ith-lead and the barrier M. 
S^'s are the self-energies that couple the scattering re- 
gion to the remaining parts of the semi-infinite leads. 
The Hamiltonian and the charge density matrix are as- 
sumed to be converged to the bulk values in the leads 
outside the scattering region. For practical calculations, 
this assumption is tested by including larger fractions of 
the leads in the scattering region and by examining the 
charge convergence during the SC iterations. Here the 
term "charge convergence" refers to the conservation of 
total charge in the scattering region. The charge conver- 
gence criterion is of crucial importance, since a failure 
to fulfill it would indicate bad numerical convergence or 
issues with the setup of the size of the extended region. 
By achieving the charge convergence, we make sure that 
all the screening takes place within the scattering region. 

The Hamiltonian matrix in Eq.[T]can be written explic- 
itly in a tri-diagonal form, as explained below. Consider 
the subregions C3 , C4 and the region M as a single region 
C, schematically shown in Fig. [2l One can rewrite the 
above Hamiltonian matrix as: 



(2) 




where 



Hr = 



Vmc 3 





We proceed with the above tri-diagonal Hamiltonian 
for multitcrminal calculations in the same way as in a 
two-terminal case^ The most time-consuming part is the 
calculation of the Green functions, i.e., the inversion of 
a matrix (eS — H), where H is the Hamiltonian matrix 
in Eq. [T]and S is the overlap matrix. For a large system, 
the matrix can be further reduced to a tri-diagonal ma- 
trix with smaller blocks. Its inversion can be done by an 



iterative method. It can also be efficiently carried out us- 
ing sparse algebra. In fact, the numerical effort to invert 
the matrix is independent of the number of terminals. 
The matrix size (or the system size) and its sparsity play 
an important role in determining the computational cost. 
B. Density Matrix 

In this subsection, we first outline the procedure 
adopted for computing the density matrix for a two- 
terminal system^ ' 10 ! 11 and then generalize it for a multi- 
terminal system. 

As explained in detail in Ref. for a two-terminal sys- 
tem, one may write the density matrix as 



D vu , = 



de 



P2j 



(3) 



P^(e) = ^[G(e)r Li (e)Gt(e)] w ,,i=l,2 (4) 

where v and u' are the indexes of localized orbitals in the 
extended scattering region, G the Green function, and 
IY^e) = i [El ; (c) — /2 is the coupling function 

for the ith-lead. = [V <7*(e) V^] is the self-energy 

of the ith-lead that couples it to the extended-scattering 
region. 

The density matrix given in Eq. [3] is general, i.e., it 
is valid for both equilibrium or noncquilibrium electron 
transport. It can be separated into two parts 
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EB 



de G(e + iS) np(e — /.ti) 



+ / dep u l,(e) (n F (e - /i 2 ) - n F (e - /ii)) (5) 



where EB is the low energy bound for the valence band. 
The first and second parts contain the equilibrium and 
noncquilibrium density matrices, respectively. 



C. Generalized density matrix for multiterminal 
junction 

Equation [5] is now generalized to the general multiter- 
minal case. The generalized density matrix is now 
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where energy EB is chosen to be low enough to include 
all of the valence bands and p m is the electro-chemical 
potential of the mth-lead. In practice, each D uu > is calcu- 
lated separately for each p m being equal to the chemical 
potential of a given lead m. All D vu i to reduce the nu- 
merical error related to the integration. 

For p m = fii , the electro-chemical potential of the lead 
i, the density matrix is 



which satisfies J2i W% VV ' = 1j with N being the number of 
leads. The weight w % vv i is chosen to minimize the numer- 
ical error in the solution^ We test the convergence by 
increasing the density of the energy mesh, thereby mak- 
ing sure that the integration yields accurate final results. 

D. Computation of the conductance 
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EB 



de G(e + iS) n F (e — pi) 



depj v ,(e) [n F {e - Pj) - n F (e - Pi)]. 



D % vv i and Aj^, arc the equilibrium and noncquilibrium 
parts of the density matrix, respectively. The integral in 
the first part of Eq. ([B]) can be carried out with complex 
contour integral technique as in Ref. S However, the 
integral in the second part, A^,,, must be calculated on 
the real energy axis with a very dense mesh. 

Because of errors related to numerical integration, the 
computed solutions of Eq. ([6]) will not produce exactly 
the same results for all i's. So, in order to minimize the 
error in the solutions, we compute the density matrix as 

a weighted sum of D vv i in the following way: 



We apply Keldysh theory for the computation of the 
conductance of the multiterminal junction. Within the 
'electron counting' picture of transport, the conductance 
G of the junction is obtained from the transmission prob- 
abilities of all scattering channels entering from one lead 
and leaving through the other^ 

G(V)=G T(V), 

where V is the applied bias voltage. The conductance 
quantum Go = e 2 /h is the inverse von-Klitzing constant 
(i.e. the quantum of resistance, Rk ~ 25.8 k£l). The to- 
tal transmittance T(V) comprises the transmission prob- 
abilities in the 'energy window of tunneling' opened by 
VM 

Once the potential profile is self-consistently deter- 
mined, the transmission spectrum from leads U to Lj 
under the external applied bias, V = pi — Pj, can be 
calculated as 



TLje,V) 



^-Tr[r Li (e)G+(e)r L .G-(e)], (8) 
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where are the advanced and retarded Green functions 
for the extended-scattering region S, and gi { is the surface 
Green function of the zth-lead. 

The current from the lead to Lj through the molec- 
ular barrier is given by 



T Uj (z,V) [/(e- Mi )-/(e-Mj)]&, 
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FIG. 3: An initial bias profile (size 60 Bohr x 60 Bohr) for 
a four-terminal junction to be applied to the system at the 
beginning of a nonequilibrium calculation. For the planar 
molecules considered here, it has been generated by solving a 
2D Laplace's equation with appropriate boundary conditions 
(see text for details). An identical bias voltage 0.8 V is applied 
through all four leads Li, l_2, L3 and L4. The leads are denoted 
by the numbers 1, 2, 3, 4, respectively, and the scale bar of 
the potential is in eV. 

where / is the Fermi-Dirac distribution. 

Although the present multiterminal NEGF approach 
was derived and carried out only within DFT, it can serve 
as a starting point for implementation of many-body 
corrections at the quasi-particl o 14 ' 15 or self-interaction 
correctio n 16 ! 17 levels. A time-dependent formulation 18 
is also possible. 

E. Finite bias 

In a two-terminal system, the initial potential for an 
applied bias can be simply a linear interpolation indepen- 
dent of the bias between the electrodes. The situation is 
not as simple in three or four terminal system. First, one 
needs to make sure that the potentials of all electrodes 
(outside the extended-scattering region S) will be unaf- 
fected by the applied bias voltage, in other words, the 
modified potential has to match at the boundary of each 
electrode and the region S. Second, the variation of the 
potential between any two electrodes through the molec- 
ular barrier has to be continuous and uniform. Third, the 
electrostatic potential in the vacuum region between two 
arbitrary electrodes has to be realistic. In order to create 
such an initial profile for the planar molecules considered 
here, we iteratively solve the 2D Laplace's equation (as- 
suming the system is in the xy-plane) in a hypothetical 
system where the extended scattering region is empty: 

d 2 V(x,y,z) d 2 V(x,y,z) _ Q 
dx 2 dy 2 

with the following boundary conditions: (a) the initial 
potential in the scattering region is zero (or may be the 
same as the potentials of the four leads), (b) the potential 
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FIG. 4: Flowchart of the self-consistent loop used to calculate 
I-V characteristics. The terminology is explained in the text. 



in every lead is unchanged, (c) the potential towards the 
vacuum region, that is, at the corners of the box, decays. 

Solving Laplace equation yields an initial bias- 
potential profile, as shown in Fig. ((3]), for a given xy-planc 
of the system. In the simple of planar molecules, we re- 
peat this image for the other planes of the 3D system. 
In more complex cases, a 3D solution of an approximate 
initial value problem would be used. We stress that the 
solution of Laplace equation is merely a starting guess of 
the bias potential that is being updated via the SC cal- 
culation. During the course of our implementation and 
testing, we found that using this solution in the first it- 
eration significantly accelerates the convergence. It is a 
very effective guess to initiating the SC iterations, as the 
solution to Laplace equation is the correct one for the 
given setup in absence of the central part. Once the po- 
tential, charge density, etc. are converged, the final result 
is independent of the initial guess. 



F. Computational Details 

The electronic properties of the tunnel junctions dis- 
cussed in Section IIIII arc obtained within the nonequilib- 
rium Green function (NEGF) approach*^ using a basis 
of optimally localized orbitals,^^ and a multi-grid ap- 
proach. The ab initio calculations for the leads and the 
molecule are performed with the O(N) method, details of 
which can be found in Ref. [||. The exchange and correla- 
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tion terms are represented in the generalized gradient ap- 
proximation (GGA).— The electron-ion interactions are 
described by nonlocal, ultrasoft pscudopotentials.— The 
surface Green functions are calculated with a transfer- 
matrix technique in an iterative schemed The potential 
and charge density in the leads are fixed to those cor- 
responding in the bulk material. The central conductor 
part includes enough "buffer layers" of the lead so that 
the potential and the charge density match at the in- 
terfaces between the conductor and leads after the SC 
calculations. The Hartree potential is obtained by solv- 
ing Poisson equation with boundary conditions matching 
the electrostatic potentials of all the leads. The gener- 
ated SC potentials and charge density serve as inputs for 
the conductance calculations. The flowchart in Fig. (j4|) 
explains the relations between the various steps in our al- 
gorithm. The computations use a massively parallel real- 
space multigrid implementation^ of density-functional 
theory DFT3£ The wave functions and localized orbitals 
are represented on a grid with spacing of 0.335 Bohr. A 
double grid technique^ is employed to evaluate the inner 
products between the nonlocal potentials and the wave 
functions, thereby substantially reducing the computa- 
tional cost and memory without loss of accuracy. 



G. Parallelization on Supercomputers 

We now describe our multi-level parallel implemen- 
tation of the multiterminal transport theory outlined 
above. First, the matrices are distributed according to 
the two-dimensional block-cyclic data layout scheme used 
by ScaLAPACK. Depending on the matrix size, one can 
use n x n processors for matrix operations (typically n 
= 1 to 4 in our applications). Second, parallelization 
proceeds over the energy points used in the integration 
to obtain the charge density matrix. Third, potentials 
and density matrices are also parallelized over the 3D 
processor grid pe x ,pe y ,pe z , where pe x x pe y x pe z is the 
total number of processors. This step of parallelization 
drastically accelerates the Poisson equation solver during 
the self-consistent iterations. Fourth, parallelization over 
the bias points is trivial and can be achieved with nearly 
100% efficiency. 

For the zero bias calculation, the computational cost 
for a multiterminal system is about the same as that for 
a two-terminal system if the number of atoms in the scat- 
tering region is the same. However, for the multitermi- 
nal system, an additional computational time is required 
for the noncquilibrium calculation. This is because, as 
the number of leads increases, the number of terms in 
the density matrix also increases (see in Eq. 6). The 
most time-consuming part of the entire computation is 
the matrix inversion needs to calculate the Green func- 
tions. This part scales nearly linearly if one takes into 
consideration the sparsity feature of the Hamiltonian and 
overlap matrices. 
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FIG. 5: (Color online) Charge convergence of the radialene 
system with SC steps at zero bias. In the inset, a schematic 
view of the central region of the four-terminal radialene junc- 
tion is shown, with the number 1, 2, 3 and 4 marking the 
positions of the leads. The size of the system is 60 Bohr x 60 
Bohr. 



III. APPLICATIONS 
A. Radialene molecule 

In order to illustrate the proposed approach for cal- 
culating the conductance of a multiterminal molecular 
device, we choose a four-terminal junction of radialene 
molecule connected to semi-infinite carbon chains as a 
first example. A schematic diagram is shown in the inset 
of Fig. O The system has C± v symmetry. Applying our 
technique to this system, we expect to see the same sym- 
metry in the converged potential profile, which should 
also be reflected in the transmission curves. 

Our noncquilibrium Green function (NEGF) 
tcchniqu o 9 ' 19 uses a basis of optimal localized orbitals^^ 
The atom-centered orbitals arc optimized variationally 
in the equilibrium geometry. In the radialene based 
system, we include 48 atoms in the calculation and 
each atom has six orbitals with the radii of 9 Bohr. 
A self-consistent calculation is carried out within an 
extended zone around the scattering region. For this 
system, the total charge is converged after 18 steps of 
the SC process, as shown in Fig. [5] The charge density 
determines the potential. 

Fig. [5] (left) shows the converged potential profile at 
zero bias 4.88 Bohr above the atomic plane. The C± v 
symmetry of the system is reflected in its potential pro- 
file. After the convergence of the charge density is 
achieved for the equilibrium density matrix, we apply 
the bias voltage through the leads. At this stage, the 
noncquilibrium part of the density matrix (sec Eq. [6|) is 
included through an iterative process. In Fig. [6] (right) 
we show the converged potential profile after applying 
an identical bias of 0.8 V through all the leads. It 
again shows the four-fold symmetry as expected. It also 
shows that after the convergence, the potential matches 
very well at the boundary of the leads and the central 
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FIG. 6: (Color online) (Left) Self-consistent converged poten- 
tial profile (same as system size, i.e., 60 Bohr x 60 Bohr) at 
zero bias of the four-terminal radialene system, plotted 4.88 
Bohr above the atomic plane. (Right) The converged poten- 
tial profile after applying an identical bias at 0.8 V through 
all the leads. Both the images are shown in same color scale 
(in eV), to compare the relative heights of the potentials. 




FIG. 7: (Color online) (Left) Converged potential profile of 
the radialene system with a non-uniform bias voltage. (The 
color scale is in eV.) An identical bias of 0.8 V is applied 
through three of the leads (denoted by 1, 2, 3) and no bias is 
applied through the fourth lead. (Right) Potential drop along 
the central line between the leads 3 and 4. 



molecule. 

We now examine the potential drop when different bi- 
ases through the leads are applied. Here, a bias of 0.8 V 
is applied at the leads Li,l_2,l-3, while the fourth lead 
l_4 is at zero bias. After solving the Poisson equation, 
we observe a uniform potential drop between the leads 
L 3 and L 4 (see Fig. (left)). For testing purposes, we 
have also plotted the Hartree potential Vh along the line 
connecting leads L3 and L 4 in Fig. [7] (right). 

Once the potential profile is self-consistently deter- 
mined, the transmission spectrum under the applied bias 
V is calculated using the Eq.[8] The transmission curves, 
shown in Fig. [H are computed for several bias voltages, 
with the voltages being the same at l_i, l_2 and L3, while 
L4 is being kept at V = 0. The left and right panels in 
Fig. [8] show the transmission L34 and L14, respectively. 
The carbon-atoms in the lead are fixed at equidistant 
bond length. It follows that the lead has metallic char- 
acter and therefore no gap appears in the transmission 
curve (as would be the case if the chain had been sub- 
jected to a Jahn- Teller distortion). With zero bias, the 
transmission curve around the Fermi level is almost con- 
stant. This nature of transmission is expected because 
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FIG. 8: (Color online) Transmission curves of the four- 
terminal radialene system with different bias voltages. The 
Ef is the Fermi energy of the lead 4. The bias geometry is 
the same as in Fig. [7] The left and right panels show trans- 
mission through the leads L34 and L14, respectively. 



the free-electron-like sp-states of carbon contribute to the 
transmission. However, with increasing bias, the trans- 
mission curve starts to oscillate. The origin of the oscilla- 
tion is related to the choice of electrode. Here we used a 
simple, very idealized carbon chain made up of 8 atoms 
per lead. In order to examine the charge convergence 
in the scattering region, one needs to define a potential 
box around the molecule, where the Poisson equation is 
solved. However, because of the small number of carbon 
atoms in the lead, the potential box needs to be quite 
large, encompassing major parts of the leads. This cre- 
ates finite-size effects, which are reflected in the transmis- 
sions showing oscillations. As the bias is increased, the 
finite size effects also increase and this is why there are 
more oscillations in the transmission curves with higher 
bias. Therefore, the these oscillations are an artifact of 
the small number of atoms in our "test" leads. However, 
realistic molecular systems, with either thicker nanowircs 
or bulk surfaces, do not show these artifacts. Our inves- 
tigations of two-terminal systems 7 ^ 27 ' 28 further demon- 
strate this claim. In addition, our ongoing investigations 
on more realistic four-terminal molecular junctions (to be 
communicated soon) are free of such oscillations. 

Note that because of the Ca v symmetry of the system, 
the transmissions L12 and L34 at zero bias are identical 
(not shown in the figure). For the same reason, the trans- 
mission contributions L13, L14, L23 and L24 are equivalent 
(not shown in figure). 

We have also computed the I-V curves for this system, 
see Fig. [5] The current contributions through leads L34 
(red line) and L14 (blue line) are obtained in the voltage 
window ±2 V . Both curves are increasing almost linearly 
because of the constant transmission around the Fermi 
energy. 
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FIG. 9: (Color online) Current-voltage characteristics of the 
four-terminal radialene junction. The bias geometry is as 
shown in Fig. [7] The current contributions through the leads 
L34, L14 are displayed. 




FIG. 10: (Color online) A schematic diagram of the crossed- 
carbon-chains system. Three cases are considered, with dis- 
tances between the chains of 7.5, 5.0 and 2.5 Bohr. 



B. Crossed carbon chains 

As a second example, we have chosen a four-terminal 
system consisting of two crossed carbon chains. Our ob- 
jective is to vertically bring the carbon chains closer and 
closer, and to see how the current varies with the distance 
between the chains, possibly leading to a crossover from 
low to high coupling. In our study, we consider three dis- 
tances between the chains, d = 7.5, 5.0 and 2.5 Bohr. To 
build the system, we include a total of 66 carbon atoms, 
with 33 atoms in each chain. Every atom has 6 basis or- 
bitals with the radius of 9 Bohr, as in previous example. 
A schematic diagram of the system is shown in Fig. [TU1 

After the charge convergence is achieved, we apply 
a 0.5 V bias through leads Li, L2, L3 and a — 0.5 V bias 
through the fourth lead L4. The converged potential and 
the charge density of the equilibrium system (i. e., at zero 
bias) is used as the initial guess for the convergence of the 
new nonequilibrium system. Fig. Qj] (left-top and bot- 
tom) shows the converged potential profiles of the system 
when the distances between the two carbon chains are 7.5 
and 5.0 Bohr, respectively. The plotting plane is paral- 
lel to the chains and passes through one of them. Both 
figures are symmetric about the line connecting leads L3 
and I-4 . To observe the potential drop along this line 
more clearly, we plot the Hartree potentials as shown in 
Fig. [TT] [(b) and (d)]. If the two chains are sufficiently 
away from each other, e.g., separated by 7.5 Bohr, the 
potential drop between two leads is smooth, and thus 
electron tunneling along a given chain will be the largest 
in this case. As we bring the chains closer, e.g., to a 
distance of 5.0 Bohr, the probability for an electron to 
tunnel from one chain to the other increases significantly. 

Fig. Q2] shows the transmission curves in all three cases 
computed through leads L34 (left panel) and L14 (right 
panel) with zero and non-zero biases. We applied the 
same bias as in Fig. 1111 The transmission of a system 
consisting of a single ideal lead must be equal to the 
number of scattering channels at the energy 



the carbon chains are far away from each another, e.g., 
at a distance of 7.5 Bohr, the overlap integral between 
orbitals on atoms situated in two different leads is close 
to zero. Therefore, each lead in the system will tunnel 
current as an isolated electrode. This is why we observe 
a constant transmission through L34, while the transmis- 
sion through L14 is almost zero, as expected. In a car- 
bon chain system, a two-fold degenerate band crosses the 
Fermi level. Consequently, there are two scattering chan- 
nels and G = 2. As we decrease the distance between the 
chains, to 5.0 and 2.5 Bohr, the orbital overlaps between 
the chains become larger and larger. This results in a de- 
crease in transmission L34 and increase for L14. The total 
transmission that would occur only through L34 in the 
bare lead case, is now partially distributed to the other 
channels. At a finite bias, the transmission curves start 
to oscillate. As explained in the Subsection MI Al these 
oscillations are due to electrons that are attracted to the 
central region, thereby conserving the static charge and 
creating an "electron-in-a-box" effect with corresponding 
standing- wave- like oscillations. 

Fig. [T3l shows the flow of current through the channels 
L34 (left panel) and L14 (right panel). We apply bias V/2 
through the leads Li,L.2,l-3 and bias — V/2 through L 4 , 
and compute the current for distances 7.5, 5.0 and 2.5 
Bohr between the carbon chains. In all cases, as the bias 
increases, the current increases almost linearly. However, 
the nonlinearity is expected to be larger in a semicon- 
ducting system. It should be noticed that the current 
contribution through the L34 channel increases with an 
increase of the distance between the chains, while the 
current flow through the L14 channel decreases. Note 
that the currents through the L14 and L24 channels arc 
identical because of the symmetry of the system. It is 
interesting to see how the current flow through a given 
channel, say L14, decays while moving the chains away 
from each other. We have fitted the current contribu- 
tions at a bias of 1.2V as 1(d) = Io e~@ d and found the 
exponential decay constant (3 to be 0.97 Bohr -1 . 
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FIG. 11: (Color online) (a) and (c) Comparison of the con- 
verged potential profiles of the crossed-carbon-chains system 
when the distances between the two carbon chains are 7.5 and 
5.0 Bohr. In both the cases, we apply same bias 0.5 V through 
the leads Li,l_2, and L3 and bias —0.5 V through the fourth 
lead L4. The plotting passing through one of the chains, (b) 
and (d) Potential drops between the leads L3 to L4 for (a) and 




FIG. 12: (Color online) Comparison of transmission curves in 
all three cases, where the distances between the chains are 7.5, 
5.0 and 2.5 Bohr, of the crossed-carbon-chain system with zero 
and non-zero biases (as in Fig. 1 1 1 p . The left and right panels 
show transmission through leads L34 and L14, respectively. 
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FIG. 13: (Color online) Comparison of I-V curves with vary- 
ing distances between the chains of the crossed-carbon-chain 
system. The left and right panels show the current contri- 
butions through leads L34 and L14, respectively, with the bias 
voltage. 



scaling method (N is the number of electrons) for comput- 
ing the electronic properties of the lead and the central 
region. This basis is used to expand the Green functions, 
the transmission function, and the charge density under 
bias, which are sclf-consistently determined via contour 
integration. The methodology is developed to scale well 
on massively-parallel computers, and should therefore be 
applicable to systems of realistic sizes. 

To demonstrate the suitability of the new technique for 
studying electron transport in multitcrminal junctions, 
we have chosen two very simple four-terminal systems 
as test applications. In the first example, a radialcne 
system having Ci v symmetry is used to test the conser- 
vation of symmetry and the numerical robustness of our 
implementation. In the second example, we have exam- 
ined the conductance properties of two crossed carbon 
chains. The I-V characteristics of the chains show the ex- 
pected trends with the changing strength of interactions. 
These demonstrations establish the general applicability 
of the method. Since our code is efficient and highly par- 
allel, we are able to deal with rather large systems. One 
such application (a four-terminal system consisting of 
an organic molecule [9,10-Bis((2'-para-mercaptophenyl)- 
cthinyl)-anthracene] connected to four gold nanowires) 
will be published elsewhere^!. 



IV. SUMMARY AND CONCLUSIONS 
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